Basic scMultiome data analysis

Vignette for integrated Seurat/Signac analysis: https://stuartlab.org/signac/articles/pbmc_multiomic.html

Overview

1. Seurat objects
   a) Create RNA assays
   b) Create ATAC assays
   c) Cell statistics
   
2. Filtering
   a) Cell filtering by HQ cells from single omics
   b) Sample statistics
   c) Cell statistics
   
3. Preprocessing
   a) ATAC
   b) RNA
   
4. Dimensionality reduction and single-cell embedding
   a) ATAC
   b) RNA
   c) Neighbor graph
   d) Non-linear dimension reduction
   
5. Save results
   
   

Input data

1. 10x Genomics Multiome-scRNA/ATAC-seq of PBMCs
   PBMC_Multiome_10xTn5_CellRangerARC_fragments.tsv.gz
PBMC_Multiome_10xTn5_CellRangerARC/raw_feature_bc_matrix
   
2. 10x Genomics Multiome-scRNA/TurboATAC-seq (1:6 concentration of Tn5) of PBMCs
   PBMC_Multiome_Tn5hc_50perc_CellRangerARC_fragments.tsv.gz
PBMC_Multiome_Tn5hc_50perc_CellRangerARC/raw_feature_bc_matrix
   
3. 10x Genomics Multiome-scRNA/TurboATAC-seq (1:3 concentration of Tn5) of PBMCs
   PBMC_Multiome_Tn5hc_CellRangerARC_fragments.tsv.gz
PBMC_Multiome_Tn5hc_CellRangerARC/raw_feature_bc_matrix
   
4. HQ cells from ATAC
   RawData_scATAC_*.csv
   For prep see 20230912_PBMC_Multiome_scATAC_UMAP.ipynb
   
5. HQ cells from RNA
   RawData_scRNA_*.csv
   For prep see 20230912_PBMC_Multiome_scRNA_UMAP.ipynb
   

Notebook

20230912_PBMC_Multiome_Co-embedding.ipynb

---

0. Load libraries, scripts and data

### Libraries
library(magrittr)
library(qs)
library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v86)
library(BSgenome.Hsapiens.UCSC.hg38)
library(ggplot2)
library(ggpubr)

### Set seed
set.seed(13)

library(future)
plan("multicore", workers = 20)
options(future.globals.maxSize = 500000 * 1024^2) # for 500 Gb RAM

### Samples
samples <- c("PBMC_scATAC", "PBMC_scTurboATAC_16", "PBMC_scTurboATAC_13")

### Sample palette
sample_palette <- c("grey", "blue", "darkblue")
names(sample_palette) <- samples

### Get gene annotations for hg38
annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotation) <- "UCSC"

### Data
inputFiles <- c("PBMC_Multiome_10xTn5_CellRangerARC_fragments.tsv.gz",
                "PBMC_Multiome_Tn5hc_50perc_CellRangerARC_fragments.tsv.gz",
                "PBMC_Multiome_Tn5hc_CellRangerARC_fragments.tsv.gz")
names(inputFiles) <- samples

### Filtered ATAC cells
atac_cells <- list()

for (sample in samples){
    atac_cells[[sample]] <- read.table(paste0("RawData_scATAC_", 
                                              sample, ".csv"), sep = ",", header = TRUE)
    colnames(atac_cells[[sample]])[1] <- "CellID_ArchR"
    atac_cells[[sample]]$CellID <- gsub(paste0(sample, "#"), "", atac_cells[[sample]]$CellID_ArchR)
}

### Filtered RNA cells
rna_cells <- list()

for (sample in samples){
    rna_cells[[sample]] <- read.table(paste0("RawData_scRNA_", 
                                              sample, ".csv"), sep = ",", header = TRUE)
    colnames(rna_cells[[sample]])[1] <- "CellID_Seurat"
    rna_cells[[sample]]$CellID <- gsub(paste0(sample, "_"), "", rna_cells[[sample]]$CellID_Seurat)
}

---

1. Seurat objects

overview

a) Create RNA assays

input_data <- list()
for (sample in samples){
    input_data[[sample]] <- Read10X(data.dir = paste0(sample, "/raw_feature_bc_matrix"))
}

proj_list <- list()
for (sample in samples){
    proj_list[[sample]] <- CreateSeuratObject(counts = input_data[[sample]][["Gene Expression"]], project = sample, assay = "RNA")
    proj_list[[sample]]$orig.ident <- sample
}

for (sample in samples){
    proj_list[[sample]]$percent.mt <- PercentageFeatureSet(proj_list[[sample]], pattern = "^MT-")
}

b) Create ATAC assays

frag_list <- list()

for (sample in samples){
    frag_list[[sample]] <- CreateFragmentObject(inputFiles[sample])
}

atac_list <- list()

for (sample in samples){
    proj_list[[sample]][["ATAC"]] <- CreateChromatinAssay(counts = input_data[[sample]][["Peaks"]], fragments = frag_list[[sample]],
                                                          genome = 'hg38', sep = c(":", "-"), annotation = annotation)
}

for (sample in samples){
    DefaultAssay(proj_list[[sample]]) <- "ATAC"
    
    proj_list[[sample]] <- TSSEnrichment(proj_list[[sample]])
    proj_list[[sample]] <- NucleosomeSignal(proj_list[[sample]])
}

proj_list

$PBMC_scATAC
An object of class Seurat 
156885 features across 720228 samples within 2 assays 
Active assay: ATAC (120284 features, 0 variable features)
 1 other assay present: RNA

$PBMC_scTurboATAC_16
An object of class Seurat 
149742 features across 716191 samples within 2 assays 
Active assay: ATAC (113141 features, 0 variable features)
 1 other assay present: RNA

$PBMC_scTurboATAC_13
An object of class Seurat 
137725 features across 719588 samples within 2 assays 
Active assay: ATAC (101124 features, 0 variable features)
 1 other assay present: RNA

c) Cell statistics

plot_list <- list()

for (sample in samples){
    plot_list[[paste0(sample, "_nCountRNA_nCountATAC")]] <- FeatureScatter(proj_list[[sample]], feature1 = "nCount_RNA", feature2 = "nCount_ATAC", 
                                                                           cols = sample_palette, group.by = "orig.ident") + 
                                                                theme(legend.position="none") + 
                                                                scale_x_log10(labels = scales::comma) + scale_y_log10(labels = scales::comma) +
                                                                ggtitle(paste(sample, " (", length(proj_list[[sample]]$orig.ident), "cells)"))
    plot_list[[paste0(sample, "_percentMT_TSSenrichment")]] <- FeatureScatter(proj_list[[sample]], feature1 = "percent.mt", feature2 = "TSS.enrichment",  
                                                                              cols = sample_palette, group.by = "orig.ident") + 
                                                                theme(legend.position="none") + 
                                                                ggtitle(paste(sample, " (", length(proj_list[[sample]]$orig.ident), "cells)"))
    plot_list[[paste0(sample, "_percentMT_NucleosomeRatio")]] <- FeatureScatter(proj_list[[sample]], feature1 = "percent.mt", feature2 = "nucleosome_signal",  
                                                                              cols = sample_palette, group.by = "orig.ident") + 
                                                                theme(legend.position="none") + 
                                                                ggtitle(paste(sample, " (", length(proj_list[[sample]]$orig.ident), "cells)"))
}

options(repr.plot.width=18, repr.plot.height=5*length(samples))
ggarrange(plotlist=plot_list, ncol=3, nrow=length(samples))

Warning message:
"Transformation introduced infinite values in continuous x-axis"
Warning message:
"Transformation introduced infinite values in continuous y-axis"
Warning message:
"Removed 149665 rows containing missing values (`geom_scattermore()`)."
Warning message:
"Removed 219953 rows containing missing values (`geom_scattermore()`)."
Warning message:
"Transformation introduced infinite values in continuous x-axis"
Warning message:
"Transformation introduced infinite values in continuous y-axis"
Warning message:
"Removed 236481 rows containing missing values (`geom_scattermore()`)."
Warning message:
"Removed 291334 rows containing missing values (`geom_scattermore()`)."
Warning message:
"Transformation introduced infinite values in continuous x-axis"
Warning message:
"Transformation introduced infinite values in continuous y-axis"
Warning message:
"Removed 247437 rows containing missing values (`geom_scattermore()`)."
Warning message:
"Removed 291442 rows containing missing values (`geom_scattermore()`)."

summary <- lapply(samples, 
                       function(sample){c(sample, '01_raw',
                                     nrow(proj_list[[sample]]@meta.data),
                                     median(proj_list[[sample]]$nFeature_RNA, na.rm = TRUE) %>% round(., 0),
                                     median(proj_list[[sample]]$nCount_RNA, na.rm = TRUE) %>% round(., 0),
                                     median(proj_list[[sample]]$percent.mt, na.rm = TRUE) %>% round(., 2),
                                     median(proj_list[[sample]]$nFeature_ATAC, na.rm = TRUE) %>% round(., 0),
                                     median(proj_list[[sample]]$nCount_ATAC, na.rm = TRUE) %>% round(., 0),
                                     median(proj_list[[sample]]$TSS.enrichment, na.rm = TRUE) %>% round(., 2),
                                     median(proj_list[[sample]]$nucleosome_signal, na.rm = TRUE) %>% round(., 2))}) %>% do.call(rbind, .)
colnames(summary) <- c('sample', 'step', 'nCells', 'nFeature_RNA', 'nCount_RNA', 'percent.mt', 'nFeature_ATAC', 'nCount_ATAC', 'TSSenrichment', 'NucleosomeSignal')

summary

A matrix: 3 × 10 of type chr

sample	step	nCells	nFeature_RNA	nCount_RNA	percent.mt	nFeature_ATAC	nCount_ATAC	TSSenrichment	NucleosomeSignal
PBMC_scATAC	01_raw	720228	2	2	22.22	2	4	0.04	1.29
PBMC_scTurboATAC_16	01_raw	716191	1	1	20.69	1	2	0	1.75
PBMC_scTurboATAC_13	01_raw	719588	1	1	12.86	1	2	0	2

---

2. Filtering

overview

HQ_cells <- data.frame(Sample = samples,
                       HQ_atac_cells = sapply(atac_cells, nrow),
                       HQ_rna_cells = sapply(rna_cells, nrow),
                       HQ_cells = sapply(samples, function(sample){sum(atac_cells[[sample]]$CellID %in% rna_cells[[sample]]$CellID)}))
HQ_cells

A data.frame: 3 × 4

	Sample	HQ_atac_cells	HQ_rna_cells	HQ_cells
	<chr>	<int>	<int>	<int>
PBMC_scATAC	PBMC_scATAC	4141	4277	3620
PBMC_scTurboATAC_16	PBMC_scTurboATAC_16	6785	7389	5999
PBMC_scTurboATAC_13	PBMC_scTurboATAC_13	5765	6529	5077

a) Cell filtering by HQ cells from single omics

for (sample in samples){
    proj_list[[sample]]$HQ_cell <- rownames(proj_list[[sample]]@meta.data)  %in% atac_cells[[sample]]$CellID[which(atac_cells[[sample]]$CellID %in% rna_cells[[sample]]$CellID)]
    proj_list[[sample]] <- subset(proj_list[[sample]], subset = HQ_cell)
}

summary <- lapply(samples, 
                       function(sample){c(sample, '02_HQcells',
                                     nrow(proj_list[[sample]]@meta.data),
                                     median(proj_list[[sample]]$nFeature_RNA) %>% round(., 0),
                                     median(proj_list[[sample]]$nCount_RNA) %>% round(., 0),
                                     median(proj_list[[sample]]$percent.mt) %>% round(., 2),
                                     median(proj_list[[sample]]$nFeature_ATAC) %>% round(., 0),
                                     median(proj_list[[sample]]$nCount_ATAC) %>% round(., 0),
                                     median(proj_list[[sample]]$TSS.enrichment) %>% round(., 2),
                                     median(proj_list[[sample]]$nucleosome_signal) %>% round(., 2))}) %>% do.call(rbind, .) %>% rbind(summary, .)
colnames(summary) <- c('sample', 'step', 'nCells', 'nFeature_RNA', 'nCount_RNA', 'percent.mt', 'nFeature_ATAC', 'nCount_ATAC', 'TSSenrichment', 'NucleosomeSignal')

summary

A matrix: 6 × 10 of type chr

sample	step	nCells	nFeature_RNA	nCount_RNA	percent.mt	nFeature_ATAC	nCount_ATAC	TSSenrichment	NucleosomeSignal
PBMC_scATAC	01_raw	720228	2	2	22.22	2	4	0.04	1.29
PBMC_scTurboATAC_16	01_raw	716191	1	1	20.69	1	2	0	1.75
PBMC_scTurboATAC_13	01_raw	719588	1	1	12.86	1	2	0	2
PBMC_scATAC	02_HQcells	3620	1844	4024	7.98	15708	43930	4.45	1.09
PBMC_scTurboATAC_16	02_HQcells	5999	1293	2253	5.64	11426	28155	4.2	1.18
PBMC_scTurboATAC_13	02_HQcells	5077	1392	2501	4.71	11431	27180	3.62	1.46

b) Sample statistics

plot_list <- list()

for (feature in c("nFeature_ATAC", "nCount_ATAC", "TSS.enrichment", "nucleosome_signal", "nFeature_RNA", "nCount_RNA", "percent.mt")){
    plot_list[[feature]] <- ggplot(lapply(proj_list, function(x){data.frame(Sample = factor(x$orig.ident, levels = samples), 
                                                                            y = x@meta.data[, feature])}) %>% do.call(rbind, .), 
                                   aes(x = Sample, y = y)) + 
                            geom_violin(aes(fill = Sample), alpha = 0.8) + scale_fill_manual(values = sample_palette) +
                            geom_boxplot(alpha = 0.001) + 
                            xlab("") + ylab(feature) +
                            theme_minimal() + theme(legend.position = "none", plot.title = element_text(size=14), 
                                                    axis.title = element_text(size=12), axis.text = element_text(size=12), 
                                                    axis.text.x=element_text(angle=45, vjust=1, hjust=1), 
                                                    legend.title = element_text(size=12), legend.text = element_text(size=12))
}

plot_list$nCount_ATAC <- plot_list$nCount_ATAC + scale_y_log10(labels = scales::comma) + ylab("log10(nCount_ATAC)")
plot_list$nCount_RNA <- plot_list$nCount_RNA + scale_y_log10(labels = scales::comma) + ylab("log10(nCount_RNA)")

options(repr.plot.width=15, repr.plot.height=12)
ggarrange(plotlist = plot_list, ncol = 4, nrow = 2)

c) Cell statistics

plot_list <- list()

for (sample in samples){
    plot_list[[paste0(sample, "_nCountRNA_nCountATAC")]] <- FeatureScatter(proj_list[[sample]], feature1 = "nCount_RNA", feature2 = "nCount_ATAC", 
                                                                           cols = sample_palette, group.by = "orig.ident") + 
                                                                theme(legend.position="none") + 
                                                                scale_x_log10(labels = scales::comma) + scale_y_log10(labels = scales::comma) +
                                                                ggtitle(paste(sample, " (", length(proj_list[[sample]]$orig.ident), "cells)"))
    plot_list[[paste0(sample, "_percentMT_TSSenrichment")]] <- FeatureScatter(proj_list[[sample]], feature1 = "percent.mt", feature2 = "TSS.enrichment",  
                                                                              cols = sample_palette, group.by = "orig.ident") + 
                                                                theme(legend.position="none") + 
                                                                ggtitle(paste(sample, " (", length(proj_list[[sample]]$orig.ident), "cells)"))
    plot_list[[paste0(sample, "_percentMT_NucleosomeRatio")]] <- FeatureScatter(proj_list[[sample]], feature1 = "percent.mt", feature2 = "nucleosome_signal",  
                                                                              cols = sample_palette, group.by = "orig.ident") + 
                                                                theme(legend.position="none") + 
                                                                ggtitle(paste(sample, " (", length(proj_list[[sample]]$orig.ident), "cells)"))
}

options(repr.plot.width=18, repr.plot.height=5*length(samples))
ggarrange(plotlist=plot_list, ncol=3, nrow=length(samples))

---

3. Preprocessing

overview

a) ATAC

for (sample in samples){
    DefaultAssay(proj_list[[sample]]) <- "ATAC"
    
    proj_list[[sample]] <- FindTopFeatures(proj_list[[sample]], min.cutoff = 5)
    proj_list[[sample]] <- RunTFIDF(proj_list[[sample]])
}

b) RNA

for (sample in samples){
    DefaultAssay(proj_list[[sample]]) <- "RNA"
    
    proj_list[[sample]] <- SCTransform(proj_list[[sample]])
}

---

4. Dimensionality reduction and single-cell embedding

overview

a) ATAC

for (sample in samples){
    DefaultAssay(proj_list[[sample]]) <- "ATAC"
    
    proj_list[[sample]] <- RunSVD(proj_list[[sample]])
}

plot_list <- list()

for (sample in samples){
    plot_list[[sample]] <- ElbowPlot(proj_list[[sample]], reduction = "lsi", ndims = 50)
}

options(repr.plot.width=6*length(samples), repr.plot.height=5)
ggarrange(plotlist = plot_list, ncol = length(samples), nrow = 1)

LSIs_list <- list(25, 25, 25)
names(LSIs_list) <- samples

b) RNA

for (sample in samples){
    DefaultAssay(proj_list[[sample]]) <- "SCT"
    
    proj_list[[sample]] <- RunPCA(proj_list[[sample]])
}

plot_list <- list()

for (sample in samples){
    plot_list[[sample]] <- ElbowPlot(proj_list[[sample]], reduction = "pca", ndims = 50)
}

options(repr.plot.width=6*length(samples), repr.plot.height=5)
ggarrange(plotlist = plot_list, ncol = length(samples), nrow = 1)

PCs_list <- list(30, 30, 30)
names(PCs_list) <- samples

c) Neighbor graph

for (sample in samples){
    proj_list[[sample]] <- FindMultiModalNeighbors(proj_list[[sample]], reduction.list = list("pca", "lsi"), 
                                                   dims.list = list(1:PCs_list[[sample]], 2:LSIs_list[[sample]]), 
                                                   modality.weight.name = "RNA.weight")
    
    for (resolution in c(0.2, 0.5, 0.8)){
        proj_list[[sample]] <- FindClusters(proj_list[[sample]], resolution = resolution, graph.name = "wsnn")
    }
}

d) Non-linear dimension reduction

for (sample in samples){
    proj_list[[sample]] <- RunUMAP(proj_list[[sample]], nn.name = "weighted.nn", assay = "RNA")
}

plot_list <- list()

for (sample in samples){
    plot_list[[sample]] <- list()
    
    plot_list[[sample]]$clusters0.2 <- DimPlot(proj_list[[sample]], reduction = "umap", group.by = "wsnn_res.0.2")
    plot_list[[sample]]$clusters0.5 <- DimPlot(proj_list[[sample]], reduction = "umap", group.by = "wsnn_res.0.5")
    plot_list[[sample]]$clusters0.8 <- DimPlot(proj_list[[sample]], reduction = "umap", group.by = "wsnn_res.0.8")
    
    plot_list[[sample]]$nCount_RNA <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "nCount_RNA")
    plot_list[[sample]]$nFeature_RNA <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "nFeature_RNA")
    plot_list[[sample]]$percentMT <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "percent.mt")
    
    plot_list[[sample]]$nCount_ATAC <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "nCount_ATAC")
    plot_list[[sample]]$nFeature_ATAC <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "nFeature_ATAC")
    plot_list[[sample]]$NucRatio <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "nucleosome_signal")
    plot_list[[sample]]$TSSenrichment <- FeaturePlot(proj_list[[sample]], reduction = "umap", features = "TSS.enrichment")

}

options(repr.plot.width=20, repr.plot.height=24)
for (sample in samples){
    ggarrange(plotlist = plot_list[[sample]], ncol = 3, nrow = 4) %>% print(.)
}

---

5. Save results

### Save Seurat object
qsave(proj_list, 'SeuratObjects.qs')

### Save Sample statistics
write.table(summary, 'Summary.csv', sep=',', col.names = T, row.names = F)

### Save raw data
lapply(samples, function(x){all(rownames(proj_list[[sample]]@meta.data) == rownames(proj_list[[sample]]@reductions$umap@cell.embeddings)) %>% print(.);
                            df <- cbind(proj_list[[sample]]@meta.data,
                                        proj_list[[sample]]@reductions$umap@cell.embeddings)
                            write.csv(df, paste0("RawData_", x, ".csv"))})

#

Isabelle Seufert
Division of Chromatin Networks, DKFZ
12.09.2023

sessionInfo()

R version 4.2.2 (2022-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /home/bq_ilander/miniconda3/envs/RWireX_v0.2.02/lib/libopenblasp-r0.3.21.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] future_1.32.0                     ggpubr_0.4.0                     
 [3] ggplot2_3.4.2                     BSgenome.Hsapiens.UCSC.hg38_1.4.4
 [5] BSgenome_1.66.3                   rtracklayer_1.58.0               
 [7] Biostrings_2.66.0                 XVector_0.38.0                   
 [9] EnsDb.Hsapiens.v86_2.99.0         ensembldb_2.22.0                 
[11] AnnotationFilter_1.22.0           GenomicFeatures_1.50.2           
[13] AnnotationDbi_1.60.0              Biobase_2.58.0                   
[15] GenomicRanges_1.50.1              GenomeInfoDb_1.34.9              
[17] IRanges_2.32.0                    S4Vectors_0.36.0                 
[19] BiocGenerics_0.44.0               SeuratObject_4.1.3               
[21] Seurat_4.3.0                      Signac_1.9.0                     
[23] qs_0.25.4                         magrittr_2.0.3                   

loaded via a namespace (and not attached):
  [1] utf8_1.2.3                  R.utils_2.12.1             
  [3] spatstat.explore_3.1-0      reticulate_1.26            
  [5] tidyselect_1.2.0            RSQLite_2.3.1              
  [7] htmlwidgets_1.6.2           grid_4.2.2                 
  [9] BiocParallel_1.32.5         Rtsne_0.16                 
 [11] munsell_0.5.0               ragg_1.2.5                 
 [13] codetools_0.2-19            ica_1.0-3                  
 [15] pbdZMQ_0.3-9                miniUI_0.1.1.1             
 [17] withr_2.5.0                 spatstat.random_3.1-4      
 [19] colorspace_2.1-0            progressr_0.13.0           
 [21] filelock_1.0.2              knitr_1.42                 
 [23] uuid_1.1-0                  rstudioapi_0.14            
 [25] ROCR_1.0-11                 ggsignif_0.6.4             
 [27] tensor_1.5                  listenv_0.9.0              
 [29] labeling_0.4.2              MatrixGenerics_1.10.0      
 [31] repr_1.1.6                  GenomeInfoDbData_1.2.9     
 [33] polyclip_1.10-4             farver_2.1.1               
 [35] bit64_4.0.5                 parallelly_1.35.0          
 [37] vctrs_0.6.2                 generics_0.1.3             
 [39] xfun_0.39                   biovizBase_1.46.0          
 [41] BiocFileCache_2.6.0         R6_2.5.1                   
 [43] bitops_1.0-7                spatstat.utils_3.0-3       
 [45] cachem_1.0.8                DelayedArray_0.24.0        
 [47] promises_1.2.0.1            BiocIO_1.8.0               
 [49] scales_1.2.1                nnet_7.3-19                
 [51] gtable_0.3.3                Cairo_1.6-0                
 [53] globals_0.16.2              goftest_1.2-3              
 [55] rlang_1.1.1                 systemfonts_1.0.4          
 [57] RcppRoll_0.3.0              splines_4.2.2              
 [59] rstatix_0.7.2               lazyeval_0.2.2             
 [61] dichromat_2.0-0.1           checkmate_2.2.0            
 [63] broom_1.0.4                 spatstat.geom_3.2-1        
 [65] yaml_2.3.7                  reshape2_1.4.4             
 [67] abind_1.4-5                 backports_1.4.1            
 [69] httpuv_1.6.10               Hmisc_5.1-0                
 [71] tools_4.2.2                 ellipsis_0.3.2             
 [73] RColorBrewer_1.1-3          ggridges_0.5.4             
 [75] Rcpp_1.0.10                 plyr_1.8.8                 
 [77] base64enc_0.1-3             progress_1.2.2             
 [79] zlibbioc_1.44.0             purrr_1.0.1                
 [81] RCurl_1.98-1.9              prettyunits_1.1.1          
 [83] rpart_4.1.19                deldir_1.0-6               
 [85] pbapply_1.7-0               cowplot_1.1.1              
 [87] zoo_1.8-12                  SummarizedExperiment_1.28.0
 [89] ggrepel_0.9.3               cluster_2.1.4              
 [91] data.table_1.14.8           scattermore_1.0            
 [93] lmtest_0.9-40               RANN_2.6.1                 
 [95] ProtGenerics_1.30.0         fitdistrplus_1.1-11        
 [97] matrixStats_0.63.0          stringfish_0.15.7          
 [99] hms_1.1.3                   patchwork_1.1.2            
[101] mime_0.12                   evaluate_0.21              
[103] xtable_1.8-4                XML_3.99-0.14              
[105] gridExtra_2.3               compiler_4.2.2             
[107] biomaRt_2.54.0              tibble_3.2.1               
[109] KernSmooth_2.23-21          crayon_1.5.2               
[111] R.oo_1.25.0                 htmltools_0.5.5            
[113] later_1.3.1                 Formula_1.2-5              
[115] tidyr_1.3.0                 RcppParallel_5.1.6         
[117] DBI_1.1.3                   RApiSerialize_0.1.2        
[119] dbplyr_2.3.2                MASS_7.3-60                
[121] rappdirs_0.3.3              car_3.1-2                  
[123] Matrix_1.5-4                cli_3.6.1                  
[125] R.methodsS3_1.8.2           parallel_4.2.2             
[127] igraph_1.4.2                pkgconfig_2.0.3            
[129] GenomicAlignments_1.34.0    foreign_0.8-84             
[131] sp_1.6-0                    IRdisplay_1.1              
[133] plotly_4.10.1               spatstat.sparse_3.0-1      
[135] xml2_1.3.3                  VariantAnnotation_1.44.0   
[137] stringr_1.5.0               digest_0.6.31              
[139] sctransform_0.3.5           RcppAnnoy_0.0.20           
[141] spatstat.data_3.0-1         rmarkdown_2.21             
[143] leiden_0.4.3                fastmatch_1.1-3            
[145] htmlTable_2.4.1             uwot_0.1.14                
[147] restfulr_0.0.15             curl_4.3.3                 
[149] shiny_1.7.4                 Rsamtools_2.14.0           
[151] rjson_0.2.21                lifecycle_1.0.3            
[153] nlme_3.1-162                jsonlite_1.8.4             
[155] carData_3.0-5               viridisLite_0.4.2          
[157] fansi_1.0.4                 pillar_1.9.0               
[159] lattice_0.21-8              KEGGREST_1.38.0            
[161] fastmap_1.1.1               httr_1.4.6                 
[163] survival_3.5-5              glue_1.6.2                 
[165] png_0.1-8                   bit_4.0.5                  
[167] stringi_1.7.12              blob_1.2.4                 
[169] textshaping_0.3.6           memoise_2.0.1              
[171] IRkernel_1.3.1              dplyr_1.1.2                
[173] irlba_2.3.5.1               future.apply_1.10.0